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Abstract. Cubature methods, a powerful alternative to Monte Carlo due to 
Kusuoka [Adv. Math. Econ. 6, 69-83, 2004] and Lyons- Victoir [Proc. R. Soc. 
Lond. Ser. A 460, 169—198, 2004], involve the solution to numerous auxil- 
iary ordinary differential equations. With focus on the Ninomiya- Victoir algo- 
rithm [Appl. Math. Fin. 15, 107—121, 2008], which corresponds to a concrete 
level 5 cubature method, we study some parametric diffusion models motivated 
from financial applications, and exhibit structural conditions under which all 
involved ODEs can be solved explicitly and efficiently. We then enlarge the 
class of models for which this technique applies, by introducing a (model- 
dependent) variation of the Ninomiya- Victoir method. Our method remains 
easy to implement; numerical examples illustrate the savings in computation 
time. 



1. Introduction 

We deal with the common problem in quantitative finance to compute, as fast 
and accurately as possible, 

(1) E[/(Xr)]- 

Here, / : M N — > M denotes a typical payoff function and (Xi) Q<t<T is an TV- 
dimensional diffusion process, given in terms of a stochastic differential equation 
(SDE) in Stratonovich form 

(YS = ij£v}(X(s,x))odBi\ 



( Xi(t,x) \ ( x\ 



J^V^X^x))^ 



\X N (t, x)j \x N ) Vjo V N (X(s, x))dsj \ZU So VfiXis, x)) ° 6B{) 

where x = (x x , . . . , x N ) € R N and B = (B 1 , ...,B d ) is a d-dimensional standard 
Brownian motion. Whenever convenient, we shall use the compact notation 

pt d pt 

(2) X(t,x)=x+ V (X[s,x))ds + y2 Vj{X(s,x))odBi, 

Jo j=1 Jo 

or, in ltd form, 

, 4 d 

X(t,x)=x + 

/.i 

3- 

where V*(x) - V %x) + | ^Ll V j° d k V H x )- 



[ V (X( s ,x))ds + V / Vj(X(s,x))dBi, 
Jo =1 Jo 
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As is common in the analysis of higher-order, weak approximation methods for 
such SDEs (cf. the classics Kloeden and Platen [7], Glasserman [4] as well as 
Kusuoka jS], Lyons and Victoir [T^] and Ninomiya and Victoir [TS] for cubature 
type methods) we shall assume that the payoff function / and all vector fields 
Vo,Vi, . . . ,Vd are smooth, with bounded derivatives of any order. The standing 
remark in this subject, implicit in all of the aforementioned references, is that any 
scheme obtained from such an analysis can and will be applied to typical financial 
diffusion models (such as Heston, SABR and their - possibly higher-dimensional - 
generalizations) even if they do not satisfy the technical assumptions initially used 
in the analysis; numerical experiments (which are necessary for every numerical 
scheme in any case!) serve as a posteriori justification^ 

We do not wish to impose any special structure on ([2]); in particular the vector 
fields are not supposed to commute (cf. Kloeden and Platen [7] [page 348] for the 
advantages in such a case in the particular case of the Milstein scheme), no affine 
structure (as in the Heston model) is assumed, nor do we want to rely on heat-kernel 
based expansions of ([l} (such as the SABR formula). In this generality, one has 
essentially two approaches. The PDE method, based on the Feynman-Kac formula, 
consists in solving the Cauchy problem for the partial differential equation 

dtu (t, x) + Lu(t, x) = 0, u (T, x) = f(x) 

where the 2nd order differential operator L is given in Hormander form L = 
Vq + \ J2i=i Vi where vector-fields are identified with first order differential opera- 
tors. As is well known, that PDE approach is prohibitively slow in higher dimension; 
there are also stability issues when L is not elliptic. The other approach is the prob- 
abilistic "simulation" method which requires two steps. In step 1 one discretizes 
X (t, x) in order to obtain an approximation X (t,x); typically, K corresponds 
to the number of partitions of [0, T]; examples include the Euler-Maruyama (EM) 
scheme 

X iEM) ' K (0,x) = 

x^' K (±±l,x) = a (bm) ^(A,^ + ^ o( ^ m) ^(A,,))x| 

where ( zf) is a family of independent Af(0, 1)0 random variables, as well as higher 

order (Milstein, Kusuoka, Ninomiya- Victoir, . . . ) schemes which we do not wish 
to detail at this moment. The discretization error is given by 




O (T/K) for Euler-Maruyama 
O ((T/K) 2 ) for Ninomiya- Victoir 



It is possible to analyze mollified/truncated versions of CIR, Heston, SABR, ...and thus 
provide further mathematical justification. For instance, it was only recently shown in full rigor 
that the classical Euler-Maruyama scheme applied to the Heston model converges; see e.g. Mao 
and Higham \E\. Let us also mention the work of Alfonsi [T] in this context. Such considerations 
are not the purpose of the present paper. 

2 Throughout the paper J\f(fi,a 2 ) denotes the normal distribution with mean fi and variance 

a 2 . 
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In step 2 one has to integrate / \ X K (T, x)j over some domain of dimension 
D = D(K) such as0 



E 



f(x K (T,x))]=[ 



F(yx,..., Vd{K)) dyx... dy D{K) . 

Here, F denotes the dependence of / [x K (T, x)J on uniform random variables, i.e., 

F (Ux, . . . , Ud(k)) — f [X K {T,x) \ for a collection (Ui, . . . , Ud(k)) of independent 
random variables uniformly distributed on the unit interval. The right-hand-side 
is approximated by Monte Carlo (MC) or Quasi Monte Carlo (QMC), essentially 
obtained by averaging M samples of F [yi, . . . ,yo(n)) ■ These samples are random 
if created by Monte Carlo (MC) and deterministic if obtained by Quasi Monte Carlo 
(QMC). In either case, we have an integration error of the form 



MC 



(/ (X K (T, x)) , M) H - E [/ (X K (T, x) 
QMC (/ (X K (T, xj) , M) - E [/ (X K (T, x) 

The central limit theorem roughly implies that MC-intcgration error is O ( 1 j\J M 
More precisely, we have in the sense of an asymptotic equality in law, 



MC (f (x K (t, x)),M 



•Af E 



-K 



flX (t,x) 



.V 



f(x K (t,x 



/M 



so that, using V / ^ 



f X K (t,x) 



V [/ (X (t, x))\ we see that the number of sample 

points M needed to attain a given accuracy (i.e. a certain e bound for the MC- 
integration error) is roughly independent of K and the discretization algorithm. 
The situation is somewhat different for the QMC-integration error. It is known that 
there exists sequences ("sample points") such that there exists C — C (f, D (K)) 
such that for all M one has 

\D(K) 



QMC (/ (X K (T, x)) , M) (w) - E [/ (X K (T, x) 



< C 



(logM) 1 



M 



In contrast to the MC case, the number of sample points M needed by QMC to 
attain a given accuracy depends heavily on the dimension of integration D (K) and, 

possibly, on the smoothness of / (^X K ^j as a function in the points yi, . . . , Ud(k)- 
Moreover, the above error estimate is known to grossly overestimate the true error 
in many cases. 

1.1. Cubature on Wiener Space. Let us briefly put the (Kusuoka-Lyons-Victoir) 
cubature method in this context. For simplicity of notation only, we consider the 
case Vq = here. A cubature formula on Wiener space is a random variable W 
taking values in the space C 1_var ([0, 1],M ) of continuous paths of bounded variation 
with values in R d such that we have 



(3) E 



odBH ■ 


■ o dB l f J 


= E 




<tj<i 






J0<tl<-<t3<l 



The dimension D(K) will depend on the method (for instance D(K) = K X d for the Euler- 
Maruyama scheme, D(K) = K X (d + 1) for the Ninomiya— Victoir scheme). 
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for all multi-indices / = 6 with all 1 < j < m, where 

to is a fixed positive integer, the order of the cubature formula. Moreover, we 
note that since the paths of the process W are of bounded variation, the integrals 
on the right hand side of ([3]) are then understood as classical Riemann-Stieltjes 
integrals. In applications, the reference interval [0, 1] in ((3|) is typically replaced 
by some (small) interval such as [0, T/K]. (Due to Brownian scaling, however, the 
problems are equivalent; in particular, a cubature formula on [0, t] is obtained by 
a scaled version of the cubature formula on [0, 1].) In the classical paper of Lyons 
and Victoir |12) the authors actually insisted that the cubature formula is discrete 
meaning that for some positive integer k, the law of W can be written as 



E 



x,s 



where Swi is the Dirac measure on Wiener space which assign unit mass to the path 
Wi (.) , zero to every other path. (Existence and explicit knowledge of cubature 
formulas is a non-trivial problem!) The idea is now to approximate the stochastic 
differential equation ([2]) for X = X (t, x) by a family of (random) ordinary time- 
inhomogeneous differential equations, 

X(t,x;W)=x + J2 /V,(X( S ,x;W0)^Md S , 

3=1 J ° ^ 

where W now denotes a cubature formula on the interval [0, i\. A stochastic Taylor- 
expansions (e.g. chapter 18 in [2] for a discussion in the spirit of cubature) shows 



as t — !• 0. Observe that in 



that E [f (X(t,x;W))] - E [/ (X (t, x))] = O (t^^j 

the case of a discrete cubature formula E [/ (X(t, x; W))] is computed exactly (no 
integration error!) by solving k ordinary differential equations. A (big) interval 
[0,T] can be handled by dividing it into K intervals of length T/K and iterating 
this procedure but now exact computation of E [/ (X(t, x; W))} requires to solve 

k + k 2 + ■ ■ ■ + k K = O (k K ) 

ordinary differential equations. When k K becomes too big one can either perform 
a Monte Carlo simulation ( "on the cubature tree" ) or resort to recombination tech- 
niques (see Litterer and Lyons |13) for the present state of art). Let us note, how- 
ever, that in many practical applications K remains small, which helps to explain 
the numerical benefits of cubature even without recombination. 

1.2. The Ninomiya Victoir (NV) Scheme. The Ninomiya- Victoir "splitting" 
scheme, introduced in |15j . is given by 

—(NV),K 

X (0,:r) 

-^{NV),K ( [k 




e^V¥^e*^X (W) ' K (!g. )X ) if A* = -1, 
e z lVJVie& v °X iNV) > K (!§,x) if A fc = +1. 



Here e v x € R N denotes the ODE solution at unit time to y — V (y) , y (0) = x and 
the probability space carries independent random- variables (Afc), with values ±1 at 
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probability 1/2, and Af (0,1) random variables {Z 3 k ). One step in the NV scheme 
corresponds actually to a (non-discrete) cubature formula of order m = 5. To see 
this, assume Vq = for (consistent) simplicity and let (bi) denote the canonical 
basis of M. d . An Revalued random path W (lu), continuous and of bounded vari- 
ation, is then created via A (o>) e {+1. 1} and d independent Af (0, 1) realizations 
Z 1 (lo) , . . . , Z d (lo). IfA(u) = -l we take W (w) : [0, T /K) -> R d , started at say, 
to move at constant speed, first an amount Z d y/T/K in b^-direction, . . . until the 
final move Z^^jT/K in bi-direction; if A (uS) = +1 the construction is similar but 
in reversed order. When Vq ^ 0, one follows the flow of the drift vector-field for 
time T I (2K) in the first and last step of the scheme; at all intermediate steps Vq is 
followed for a time T/K; this is inspired by classical splitting methods in operator 
theory. Let us also note that the coin-flipping corresponds to Talay's trick of, in 
a weak approximation context, replacing the (difficult to sample) Levy's area by 
a discrete moment-matched random variable, see Kloeden and Platen [7] [page 466 
f.]. 

The NV scheme has attracted wide attention since its introduction in [T5] ; it is 
nowadays found in various sophisticated numerical packages such as Inria's software 
PREMIA for financial option computations. A variation of the scheme designed 
to deal with degeneracies arising some affine situations is discussed in pQ. Let us 
also mention the "NV inspired" schemes developed in 3 and [16j . 

1.3. Semi-closed form cubature. It is clear from the preceding discussion that 
cubature methods, and the NV scheme in particular, heavily rely on the ability 
to solve, fast and accurately, ordinary differential equations. The general cubature 
methods involves time-inhomogeneous ODEs; in general, there is no alternative to 
solve them numerically, typically with Runge-Kutta methods. (A detailed discus- 
sion on how Runge-Kutta methods are applied in this context is found in Ninomiya 
and Ninomiya [T4].) 

On the other hand, the Ninomiya- Victoir splitting scheme only involves the 
composition of solution flows to time-homogeneous ODEs. In particular, there will 
be "lucky" cases of models where all (or at least most) ODE flows can be solved 
exactlylj In such a case one has effectively found a level-5 cubature method which 
can be implemented without relying on numerical ODE solvers. In particular, one 
expects the cubature methods to perform especially well in such cases. As was 
observed in [TS], see also Section |2~T1 the Heston model is such a lucky case. We 
thus propose the following definition. 

Definition 1. A diffusion model of type ([2]) where a cubature method can be 
implemented without any numerical ODE solutions is said to be accessible to semi- 
closed form cubature (SCFC). 

For instance, any model of type @ where all ODE flows e tv ° , . . . , e tVd can be 
solved in closed form falls in this class. However, one soon encounters model (e.g. 
the popular SABR model, see Section [2~2j) in which some of the vector- fields do not 
allow for flows in closed form. The contribution of this paper, beyond suggesting the 

^As of Sep 2010, the weblink ralyx.inria.fr/2006/Raweb/mathfi/uid21.html contains some rel- 
evant information. 

^ By this we mean a closed-form solution to the ODE y = V (y) , y (0) = x which allows for 
fast numerical evaluation. In particular, we are not interested in "closed-form" solution in terms 
of complicated and slow-to-evaluate special functions. 
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systematic use of financial models that are accessible to semi-closed form cubature, 
is that the class of such models can be significantly enlarged by working with an 
almost trivial modification of the NV scheme. @ Before explaining our modification 
we point out that the SABR model then becomes accessible to semi-closed form 
cubature. Our modification is based on the trivial equivalence of ([2]) with 

/ d \ d 

dX(t,x)=\V (X(t,x))-^2 7j V j (X(t,x)) &t + Y J V ] {X{t,x))o&(B{+ 1] t 



^ ] (X(t, x)) + ]T V^Xit, x)) o d (b{ + ljt ) 

3 = 1 

whatever the choice of drift parameters 71, . . . , 7^. Assume that all diffusion vector- 
fields (Vi, . . . ,Vd) allow for flows in closed form, whereas e tVo is not available in 
closed form. The point is that, in a variety of concrete examples, one can pick drift 

(t) 

parameters 71 , . . . , 7<j in a way that e tv ° can be solved in closed form after all. 

Therefore, we propose the following variant of the Ninomiya-Victoir method 
(which shall be referred to as the "NV scheme with drift (trick)"): 

X iNV ">' K (0,x) = x€R N , 
(4) 

-=(NVd),K ( (k+ 1)T 
X — ,x 



(5) 



K 

e&^Vi* • • .e z t v ^ v ^X {NVd) ' K (f ,x) , if A, = -1, 
e^\ Zi v d ... e zlv le ^v^x (NVd) ' K (f ,z) , if A, = +L 



where Z\ ~ N (^7i, j?) independent of each other. 

The bulk of this paper is devoted to implement these ideas for a handful of 
(stochastic volatility) models encountered in the financial industry. Since high- 
dimensional problems are the raison d'etre for probabilistic simulation methods, 
a detailed discussion of a higher-dimensional (SABR-type) model is included. At 
last, we discuss numerical results obtained with our "drift-modified" NV scheme: 
relative to the classical NV scheme we observe significant and consistent savings in 
computational time. 

Note that we want to concentrate on the method itself, without further improve- 
ments like variance reduction, optimization of code and Romberg extrapolation. 

Acknowledgment: Partial support of MATHEON and the European Research 
Council under the European Union's Seventh Framework Programme (FP7/2007- 
2013) / ERC grant agreement nr. 258237 is gratefully acknowledged. 



^While SCFC corresponds to the "luckiest" case of avoiding numerical ODE solvers altogether, 
any significant reduction of numerical ODEs to be solved will be desirable. Our modification of 
the NV scheme can obviously be used to this purpose as well. 
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2. Application of classical NV scheme to Heston and SABR 

2.1. Heston model. The stochastic volatility model of Heston is given by the 
SDE: 



dXi(t,x) =pX x (t,x)dt+ ^JX 2 (t, x)Xi (t, x)dB] 

dX 2 (t, x) =k{6 - X 2 (t, x))dt + £y/X 2 (t,x)d [pB] + y/l - p 2 B^ , 

where p is the rate of return of the asset, 9 is the long vol, K is the mean-reversion 
rate, £ is the vol(atility) of volatility ) and p is the correlation parameter between 

the (standard) Brownian motions B\ and (pB\ + ^T~p*Bf^j . 
The vector fields are given by 



and so we get 



,vi{x))-\vi{x)) [kEUWMJ'y <o- X2 )-ie 

The corresponding solutions to the ODEs arc (cf. Lord et al. [TTJ p. 8-9] and their 
reference to [T5]; see also the Appendix) 

'xi exp ([p - \ip - \J]s + i^_2[c— - I})" 
(x 2 - J)e- KS + J 



e sVo x 



e sV >x 



Xl exp y 



Xi 



e sV *x= I r. 



with J = - — - A ^-. We assume (as in [TT]) that J > 0; see [T] for how to proceed 
otherwise. 

The Heston model can be rewritten in log-coordinates. Define Y\ (t) = log X\ (t, x) 
and Y 2 (t) — X 2 (t). In this new coordinate chart, the vector fields are 

and the corresponding solutions to the ODEs are 

e .v 0t , = (ttt + [M - ltp-U]s + §^[e— - 1 
V \ (2/2 - J)e~ KS + J 

(sips+Vyi)!.- 

e^ y =[^ + ^r— 

(2^ + 7^2)+ 
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k8 — if 

with J = — -2 — , as before and y = (2/1,2/2) = (log^i,^). As is well-known, it 
is far preferable to use Heston in log-coordinates when simulating with the EM 
scheme. Although this is less critical in the cubature context, we still recommend 
© to avoid the numerical evaluation of exp(-). 



2.2. SABR model. The SABR model is given by 

AXi{t,x) =aX 2 {t,x){Xi{t,x)) f3 dB 1 t 
dX 2 (t,x) =bX 2 (t,x)d (pB\ + y/1 - p 2 B 2 ^ , 

where h < < 1, a, b > and —1 < p < 1.0 The corresponding vector fields are 



V (x) 
and so we get 



by/l- p 2 X 2 



l^iV^Hx)] (-l[a^xlxf- l + abpx 2 x{] 



The solutions to the ODEs corresponding to the vector fields V\ and V 2 are 



e^x = 



9i(s) 
x 2 exp (bps) 



e sV2 x 



Xl 



x 2 exp I by/l — p 2 s 



where 



9i(s) 



(l-P)^(e b »°-l)+x{-e 



1 1/(1-/3) 



t \ (ax 2 , b 

gi(s) =a;iexp I — (e 



1 



+ 

13 = 1. 



< < 1, 



e sVo x 



For details on the uniqueness of g\ we refer to the Appendix. Concerning the 
solution to the ODE corresponding to Vq, let H(s) be the first component of e sVo x, 
i.e. 

H(s) 
K x 2 exp (-^b 2 s) / 

It is impossible to find H in closed- form (unless p = 0or/3 = l). This means that 
applying the standard NV-scheme must involve the numerical solution of auxiliary 
ODEs. We shall see later that with the NV scheme with drift all involved ODEs 
can be solved in closed form. 



7 Although in the literature the SABR model is also considered for < f) < i we restrict 
ourselves to the case ~ < /3 < 1 in order to avoid difficulties regarding well-posedness of X, cf. 
®. 
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3. Models accessible to SCFC and NV with drift 



3.1. Motivation. In the classical NV scheme, only centered Gaussian (Brownian) 
increments are used to flow along the diffusion vector fields. Our main observation 
is that one can also use non-centered Gaussian increments; this affects the drift 
term and, chosen in a smart way, can sometimes render all auxiliary ODE to be 
solvable in closed form. To motivate the class of models for which this works, we 
illustrate how to systematically construct models accessible to SCFC from a fairly 
general two-factor stochastic volatility model given in Ito form by 



where (Xi(0),X 2 (0)) = (x\,x 2 ) is kept fixed. In Stratonovich form this becomes 
(omitting the dependence on t in the drift and diffusion coefficients), 



In the subsequent analysis we shall exhibit a number of possible choices which lead 
to models accessible to SCFC. First we would like to choose the coefficients A,...,E 
such that we can rewrite the first SDE as 



for a constant 71 and functions Hi,H 2 . Three possible ways to achieve this goal 
are 



Note that the Heston model is a particular example satisfying (i). However, in case 
(i) we would have 71 = and since we want to illustrate the additional benefit 
of the NV scheme with drift over the classical NV scheme, we will not consider 
case (i) in any more detail. Moreover, since we would like the volatility factor to 
depend on the Brownian motion driving the stock, we will also skip case (ii) and 
concentrate on case (iii) which implies g/^) ^ -^(-^2) and 71 > 0. E.g. if we 
choose D{x) = 1, then B(x) cx exp(cx), if D(x) = x, then B(x) cx x c with c 7^ 0, if 
D{x) = x q , < q < 1, then B{x) oc exp(y^x 1_<^ ). We focus on the most natural 
choice (i.e. B not being of exponential type) and therefore pick 



dXi(t) =A{X 1 {t))B{X 2 {t))dB 1 t 

dX 2 {t) =C(X 2 (t))dt + D{X 2 {t))dBl + E{X 2 (t))dB 2 , 



dXt(t) = - \A{X X ) [A'{X 1 )B 2 (X 2 ) + D(X 2 )B'(X 2 )] dt + A{Xx)B{X 2 ) ° dB] 
dX 2 {t)= C(X 2 )-^[D(X 2 )D'(X 2 )+E(X 2 )E'(X 2 )] dt 
+ D(X 2 ) o dBl + E ( x 2) o dB 2 , 



dX^t) - F 1 (X 1 )F 2 (X 2 )di + A{X 1 )B(X 2 ) o d (Bl + 7l i) 



(i) A'(Xi) cx 1, (ii) D(X 2 ) = or (iii) D{X 2 )B'(X 2 ) cx B(X 2 ). 



B(X 2 ) = al 2 " and D(X 2 ) = b P X 2 . 

These choices give us 




With 
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define h(t;x2) = X2B tv ° ■ We would like f* h(s; X2) 2a ds to have an explicit ex- 
pression, since it will appear in the first component of xe tv ° . Possible cases are (i) 
h(t;x 2 ) oct+x 2 , (h) h(t;x 2 ) oc e ct or (hi) very specific cases like h(t; X2) = p(t, x 2 )^ 
with p nice. Both (i) and (ii) lead to C being affine and E being affine or of square 
root type. We pick 

C(X 2 ) = k(0-X 2 ) and E(X 2 ) = by/l- p 2 X 2 . 
and we shall later motivate why we let E be linear. With these choices we can write 

dXi(i) = - iA(Xi)A'pfi)a 2 Xf Q di + A(Ai)aA 2 " o d (b\ - ^abpt 



dX 2 (t) = 



k{6 - X 2 ) - X -b 2 X 2 



dt + bpX 2 o dBl + - p 2 X 2 o dBf. 



We now rewrite the second SDE in the form 



dX 2 (t) = H(X 2 )dt + +bpX 2 o d (Bl + 7l t) + by/1 - p 2 X 2 o d (B 2 t + 72 <) 

with 71 = — T^abp (as in the first SDE) and with 72 a constant such that H(X 2 ) 
becomes as simple as possible (recall that we want J*(x2e sH ) 2a ds to be explicit). 
We get 



dX 2 {t) = n6dt + bpX 2 o d ( Bl - i 



2 



abpt 



, P. ./^ abp 2 -2n/b-b 

+ by/l - p 2 X 2 o d \B 2 H ' 



Note that if we would have chosen E to be affine but not linear or have chosen E of 
square root type, we would not be able to make H so simple. (For the same reason 
we have chosen D linear and not generally affine.) 

Finally, A(Xi) is left to choose. Since we want to end up with a model accessible 
to SCFC, the function A should be such that x\e tA is explicit, which means that 
we want f x to have an explicit inverse. Also, x\e tA ' A should be explicit. The 

obvious candidates are A{X X ) = X{, A(X X ) = e cXl and A{X X ) = X x + c and all 
lead to models that are accessible to SCFC. As a case study we choose the first one 
and apply the NV scheme with drift to the resulting model in the next section. 

3.2. Generalized SABR (with shifted log-normal 2nd factor). In the pre- 
vious section we constructed a particular class of SV-models which are accessible 
to SCFC, namely: 

dX^t) =aX 2 {t) a X 1 {tfdBl 



dX 2 (t) =k{6 - X 2 (t))dt + bX 2 {t) (pdBl + v/1 - P 2 dB._ 



with ^i(O) = x\ and ^2(0) = X2- We assume that the parameters satisfy i < 
j3 < 1, 9, K > 0, a > 0, a, b > 0, — 1 < p < 1. Details surrounding well-posedness, 
intcgrability properties and martingale properties can be found in Lions and Musiela 
PHO] . A simple application of Ito's formula shows that 

X 2 (t) = x 2 e-^+i b ^ t+hWt + k9 ^ e -(-+^Kt-s) e b(w t -w s ) ds ^ 
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where W t = pB\ + y/l - p 2 B 2 and thus X 2 (t) > for all t > 0, provided x 2 > 0. 

We shall now give all the ODE solutions that are required to apply the NV 
scheme (with drift). First note that the vector fields V\ and V 2 corresponding to 
B 1 and B 2 are given by 



V x = ( "f 2 ^ ) , V 2 



bpx 2 J ' 2 \b\/T^ p*x 2 



„2 , 



and the ODE solutions are 



e sV >x 



9i(s) 
x 2 exp (bps) J ' 



/ Xi 

e sV2 x ' 



l<e<i, 



x 2 exp (byjl - p 2 s S j J ' 

where 

= [(1-^)^(6^-1)+^ 

5l ( S )^ 1 exp(^|(e— -1)), /? = !. 
The Ito drift vector field Vo and Stratonovich drift vector field Vq of X are given 

by 

V (x) = ( „ ), V (x) = (-Wftf- 1 -J-abpx^\ 
uv ; \k6-kx 2 )' uv ; ^ K 6 - (k + \b 2 )x 2 J 

We have e sVo x = (H(s),h(s)) T with 



h(s) 



V 2 K+\b 2 ) e + K+\b 2 



and H needs to be numerically solved. (We have already pointed to this difficulty 
when we discussed the classical SABR example earlier on.) 

Let us now show that by using Brownian increments with drift, this problem can 
be resolved: all necessary flows can be computed in closed form. Recall from the 
previous section that we can rewrite X as 

dX x {t) = - )^^X 2a Xf- x dt + aX^Xl o d (B\ + 7i t) 



dX 2 {t) =K0dt + b P X 2 o d (B] + 7i t) + by/1 - p 2 X 2 o d (B 2 + l2 t) . 

with 

1 , , abp 2 -2K/b-b 

7i = --abp and 72 = — . 

^ 2y/l - p 2 

We see here that the assumption — 1 < p < 1 is crucial. Note that the vector fields 

corresponding to B\ + "fit and B 2 + 7 2 i, respectively, are V\ and V 2 . Denote by 

V the remaining part, i.e. 



V^'(x) = 



k6 
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Then we have e sV ° T x = (go, n9s + x 2 ) T with (cf. Appendix) 

9o(s) = (-a 2 (3(1 P)P(s) + 4 (1 , |</?<1, 

g (s) =a?i cxp ^-i a 2 P(s)^ , f3 = 1, 
3o (s) =- ia 2 P(s) + a; 1 , /3 = i, 

where 

^)=(^Tp((^ + ^) 2a+1 -^ +1 )- 

Note that when k = or 8 = 0, P(s) should be understood in the limiting sense, 
i.e., P(s) = sx 2a for k — or 8 = 0. 

Remark 2. Since the SABR model is a special case of the model presented here 
- corresponding to a = 1, k = - the semi-closed form NV algorithm developed 
above can be, in particular, applied to the SABR model. 

3.3. Girsanov transform. We have seen for the example above, that if one uses 
the standard NV scheme, the flow of the drift vector field is not available in closed 
form. Besides using the 'drift trick' as we did above, it is also possible to absorb 
this drift in a change-of-mcasure; the details of this are outlined below. There 
is, however, a serious downside to this: the Girsanov-density which appears due 
to the change-of-measure will add significantly to the variance of the object to be 
sampled. Thus, without further variance reduction, we do not advertise the use of 
the Girsanov transform in this context. 

Let Y = (Yi,Y 2 ) be the process defined by 

dYx(t) =aY 2 (t,x) a Y 1 (tf (- 7l di + AB}) 

AY 2 (t) =k(6 - Y 2 (t))At + b P Y 2 (t) (- 7l dt + ABl) + b^l-p 2 Y 2 (t) (- 72 di + dB 2 ) . 

and let P be the probability measure under which B = (B 1 , B 2 ) is a 2-dimcnsional 
standard Brownian motion. Define the probability measure Q by ^\? t — £(t), 
where 

£(t) = cxp (jiBj + l2 B\ - i( 7l 2 + 7 2 )^ . 

Then by Girsanov, under Q, (B 1 , B 2 ) is equal in law to a 2-dimensional standard 
Brownian motion plus constant drift equal to ( 7 i, 7 2)- Hence under <Q>, Y is equal 
in law to X under P and we have for / : R 2 — > R measurable, 

E p [f(X(t))} =E®[f(Y 1 (t),Y 2 (t))} =E P lf(Y 1 (t),Y 2 (t))£(t)}. 

Hence a "weighted" NV scheme with explicit solutions to all ODEs can be ob- 
tained by using the NV scheme for the process Y and then multiplying the payoff 
f(Yi(t),Y 2 (t)), as is done in importance sampling, by £(t). Note that all the ODE 
solutions corresponding to the NV scheme for Y are explicit, since the vector fields 
corresponding to Y are Vq 7 \ V\ and V 2 . 

To back up our claim about the additional variance caused by the Girsanov 
density £(t), note that V(£(i)) = e' 7l+72 '* — 1, which is only negligible when 71 
and 72 are close to zero. 
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3.4. A multi-dimensional version. Let us illustrate how the techniques intro- 
duced (until now in the context of 2-dimensional models) remain feasible in typical 
higher dimensional models (what we have in mind here is some multi asset SV 
model) . Since it is precisely the curse of dimensionality that forces one to use sto- 
chastic methods (rather than PDE methods, say) we want to be fully explicit in 
showing how our ideas are implemented in higher dimensions. More specifically, 
we shall consider the following multi-dimensional version of our SABR-type model: 
for i = 1,...,N, 



dXi(t) =a i Y i {t) ai X i {tf i dB\ 

dYi(t) =Ki{6i - Yi{t))dt + biYi(t)dW;, 



with Xi(Q) = Xi, Yi(0) = yi and \ < fa < 1, > 0, a; > 0, a,,^ > 0, 

-1 < Pl < 1. Here (1), with B = (B 1 , . . . , B N ) T and W = (W\...,W N ) T is 
a 2N-dimensional Brownian motion with correlation matrix given by p which we 
assume to be positive-definite. Let ^fp be the unique lower-triangular matrix such 

that ^fp^fp F = p (Choleski decomposition). Then (~) ^— \/^(w)' wnere (w)> 
with B = (B 1 , . .., B N ) T and W = {W 1 , W N ) T is a 27V-dimensional standard 
Brownian motion. Hence we can write for i = 1, . . . , N, 

(N N 
3 = 1 3 = 1 

(N N 
E vpn^bi + e vPN+i,N +j wi ] ■ 
3=1 3=1 



Let for j = 1,...,N, Vj and Uj be the vector fields corresponding to B 3 and W J , 
respectively. We have 



Vj(x,y) 



( a iy ^x^^p hj \ 

b iyiVp N +i,j 

aNV% N x f) N N \fp Nj 
\ b N y Ny /p 2Nj ) 



, Uj(x,y) 



^ 

b iyiVp N +i.N+ 3 



\ b NVN^P 2NN+ J 



It follows that 



\pxa\fay\ 



2«i 201-1 



Ki0i - (ki + ^bfri)yi 



V (x,y) = 



iPNa 2 N p N y% aN x 2 N ^ 



knOn ~ (kjv + \b 2 N r N )y N 



\qNaNa N bNVN N 



J 
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where for i = 1, . . . , N, 

N r 



Pi = 



3=1 
JV 



= Pi,i = 1, 



JV 



1i [v^>JV-M,3'V^>i,j + V^jV+i,JV+jV^i,JV+j] = X/ ^PN+i,jVPi,j 

3=1 
N 

3=1 

We would like to write the system (Xi, Yi), i — 1, . . . , N in the following way 
dJQ(t) - ' 



N+i,j 



N+i,N+j 



3=1 

= PN+i,N+i = 1- 



v i=i 

dYi(t) ^Oidt + biYiit) 

JV 



' N N 

E Vp id ° d H + 7,«} + E ^p-,n +j ° d { + ¥} 



i=i 



X 



E y/Pw+v ° d H + 7;-*} + E v^jv+i,JV +J - ° d { W t + ¥} > 



3=1 



for a certain 7 = (71, • . . ,7jv) an d <5 = (Si, . . . , <5jy). Looking at Vb, we should 
choose 7 and (5 such that for i = 1, 

JV 



.,N, 

JV 



Ev^^+Evpi.jv+j^ 



1 



-AiOiih, 



3=1 



3=1 



N 



N 



E v^jv+i,j7j + E ^ 



±& 2 

2 1 



3=1 3=1 

These are 2N linear equations with 2N unknowns. It follows that there exists a 
unique 7 and 5 such that the above equalities are satisfied if yfp is of full rank. 
This is the multi-dimensional analogue of the condition — 1 < p < 1 of the previous 
section. We see that in order to apply the NV scheme with drift, we need to find 
the flows corresponding to the vector fields Vj(x,y), Uj(x,y) and 



V 



1„2 o „2a« 2/3„-l 
'2 N""i)N x jv 



All the solutions to these ODEs can be found explicitly as in Section 13.2 



3.5. Numerical analysis of our NV scheme with drift. In this section we want 

as K —> 00 

—(NVd),K 



to prove second order weak convergence of E (f (x NVd,K (T,x] 

for smooth / with A v " ' (T, x) given by (|4])-(j5]). As in the original proof by 
Ninomiya and Victoir [15] for the classical NV scheme, we use a Taylor expansion 
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to get the local order of the weak error by comparison with the known local weak 
order of the classical Ninomiya-Victoir scheme. Let X^ NV ^' K (T/K,x) be as in 



Section [l.2l Then the difference in the Taylor expansion of the expectation in one 
step is given by 



(7) 



E ( / (x (NVd) ' K (T/K,x) 



Ai = -l)-E(/fT (W (T/X,,; 



Ai = -1 



E 7i7iW(») 

l<i<j<d 



ia j v i v j f(x)+ 

i 

+ \ E 7^(^) 2 /W-J E lMV J ff{x) + 

l<i<j<d l<j<i<d 

+ 7 E TiW^/^J-i E 7^) 2 ^/(*) 



l<j<i<d 



(T/i^) 2 + O ({T/Kf\ 



\<i<j<d l<j<i<d 

When we condition on Ai = 1, only the order of the indices are swapped. By the 
peculiar structure of |0, this means that the signs of all terms of the difference 
change, i.e., 



E(f(x iNVd) ' K (T/K,x) 



Ai = l - 



■('( 



X [T/K, x) 



A^-lj-Ef/fX^V/^,,) 



Ai = 1J = 
) | Ai = -1 
+ O ({T/K) 3 ) . 



Thus, taking the unconditional expectation in Ai gives 

f— (NV),K 



E(f(x' NW (T/K,x) 



ElflX 



(T/K,x)))=o((T/Kf) 



and second order convergence of the Ninomiya-Victoir scheme with drift follows 
exactly as in the case without drift. 

Remark 3. The drift trick also works for classical cubature on Wiener space as in 
Lyons- Victoir [12], i.e., when W°(t) = t, simply by adding the same drift tji to 
the i th component of the cubature path W. This procedure retains the original 
order of convergence for the given cubature formula, as can be trivially seen by 
comparing the ODEs with and without drift. Note that in the Ninomiya-Victoir 
scheme a slightly more difficult argument as discussed above is necessary, since here 
the "time" component W° is not linear. 

4. Numerical results 

In this section we report the results of our numerical experiments. For this we 
have chosen three models: the SABR model, the generalized SABR model and 
the multi-dimensional generalized SABR model. The numerical results for these 
models are given in Section 14.11 Section 14.21 and Section 14.31 respectively. For each 
of the models, we compare the NV scheme with drift to the regular NV scheme and 
the Euler scheme and in all the experiments reported in this paper, we used Quasi 
Monte-Carlo for the integration. In order to check the order of the three schemes, we 
first give plots of the relative discretization error against the number of time steps 
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K. For these plots the number of 'simulated trajectories' M is chosen such that 
the integration error is negligible compared to the discretization error. (Since we 
are using Quasi Monte Carlo, M is strictly speaking not the number of simulations, 
but the size of the Unite low-discrepancy sequence used for the computation.) In 
order to compute these errors we of course need to know the 'true value' (or a 
good estimate of it). For the two experiments involving the generalized SABR 
model, the true value was obtained by running the code longer than reported in 
the plots. However, in the case of the SABR model, the true value was estimated 
by extrapolation of the values obtained by the code (the SABR formula is not exact 
enough!) 

Results comparing the discretization error against the number of time steps, 
might not really seem practically relevant, since they only relate the number of time- 
steps for different numerical methods, but not the corresponding computational cost 
or computer time. Therefore we also give tables of the computational time of the 
different schemes. If we want to compare the run-time for different methods, we 
need to ensure the fairness of the comparison, i.e., we need to compare the different 
methods with parameters giving similar computational errors. We basically have 
two parameters for the numerical method, namely the number of time-steps K 
and the number of simulated trajectories M - where every trajectory is given by a 
D(if )-dimensional vector from a low discrepancy sequence, in our case the Sobol- 
sequence. (Here, D{K) depends on both the model and the method.) Unlike 
for Monte Carlo simulation, where accurate, but probabilistic error estimates are 
available, there is no simple error estimation procedure for QMC (as far as we are 
aware). Therefore, it is not obvious how to choose the number of trajectories for 
a comparison of run-times. Choosing the same number of trajectories for every 
method might not be appropriate, because some methods might yield "rougher" 
integration problems, requiring a higher number of trajectories for a comparable 
precision. 

In our comparison we proceeded as follows. For a given model and method 
we first choose K (and M very large) such that the relative discretization error 
is around 10~ 3 . For simplicity, we take the computations run for producing the 
plots, which means that if is a power of two. Then wc start with Mq = 1000 
trajectories, run the method, and double the number M of trajectories until the 
observed (absolute) error is consistently closer than 2 x 10~ 5 to the true (absolute) 
error for the fixed given K. 

To explain this in more detail, first recall that the computational error in (Quasi) 
Monte Carlo simulations for SDEs splits into two parts: Error = Errors + Errorj nt . 
Here, Errors stems from the time-discretization of the SDE. This part of the error 
is controlled by K . Error; nt is the error from the integration, i.e., from the numerical 
computation of the expectation of the solution to the discretized SDE. This error 
part is controlled by M. In the comparison procedure, we first fix K such that the 
relative discretization error Errordi s /C (C denoting the true result) is around 10~ 3 . 
Then we choose M (by a doubling procedure) such that Errorj nt < 2 x 10~ 5 . (Since 
in all cases C w 0.1, this means we choose the integration error to be one fifth of 
the discretization error.) 

In the following tables, we report the found parameters, the corresponding rel- 
ative error and the computational time in seconds - all computations where per- 
formed on the same computer, a Toshiba laptop with 6 GB RAM and four Intel 
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Core i7 CPUs with 1.6 GHz. Notice that the run-time scales linearly with the 
number M of trajectories. 

4.1. SABR. In this section we give the results corresponding to the experiment 
with the SABR model. Recall that this SV-model is given by 

dXi(t) =aX 2 (t)X 1 (tfdBl 

dX 2 (t) =bX 2 (t) (p&B] + y/1 - p 2 dB f 2 ) , 

with Xi(0) = X\ and ^(0) = x 2 and where the parameters satisfy | < /3 < 1, 
a, b > 0, — 1 < p < 1. The parameters chosen for the experiment are /? = 0.9, 
a = 1.0, b = 0.4, p = —0.7, .T! = 1.0, x 2 = 0.3. As the derivative we choose a 
(European) call option with maturity time T = 1.0 and strike price K = 1.05. For 
simplicity we assume that the interest rate is zero. The corresponding estimated 
'true result' is 0.09400046. 

In Figure [T] the convergence rates of the three schemes is graphically displayed. 
We clearly see the second-order convergence of the two NV schemes compared to 
the first-order convergence of the Euler method. 




~i 1 1 1 1 1 1 1 r~ 

1 2 5 10 20 50 100 200 500 



Number of timesteps 



Figure 1. Order of convergence for the SABR model. 



Method 


K 


M 


Rel. Error 


Time 


Euler 


32 


512000 


0.00150 


5.87 sec 


Ninomiya-Victoir 


2 


512000 


0.00134 


2.44 sec 


NV with drift 


2 


128000 


0.00140 


0.28 sec 



Table 1 . Computational time for the SABR model 



In Table [T] the timings are reported for the SABR model. Notice that the 
Ninomiya-Victoir method both in its original form and in its variant are clearly 
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more efficient than the Euler method, by a factor two or three. On the other hand, 
the simpler structure of the Ninomiya-Victoir method with drift results in a con- 
siderable speed up if compared with the original method. Notice that a speed up 
with a factor two would even hold if we reject the empirically found choice of M 
for the drift variant and use the same number of trajectories as for the other two 
methods. 

4.2. Generalized SABR. In this section we give the results corresponding to the 
experiment with the generalized SABR model of Section 13.21 For convenience we 
restate this SV-model: 

dXi(t) =aX 2 {t) a X 1 {tfdB 1 t 

dX 2 (t) =k(6 - X 2 (t))dt + bX 2 (t) (pdB} + y/l - P 2 dB^ , 

with Xi(0) = xi and X 2 (0) = x 2 and \ < p < 1, 6 > 0, k > 0, a > 0, a, b > 0, 
— 1 < p < 1. For our experiment, we choose the parameters as follows: /3 = 1.0, 
9 = 0.3, K = 2.0, a = 0.5, a = 1.0, b = 0.5, p = -0.7, x x = 1.0 and x 2 = 0.2. We 
further pick the same call option as in Section fOI The estimated 'true result' is 
0.1767505855. 

In Figure [5] the discretization error against the number of time steps is plotted. 
Again, we see the second-order convergence of the two NV schemes compared to 
the first-order convergence of the Euler method. 




1 2 5 10 20 50 100 



Number of timesteps 

Figure 2 . Order of convergence for the generalized SABR model. 

In Table [5] the timings are reported for the generalized SABR model. In the 
one-dimensional generalized SABR model, the Ninomiya-Victoir method is only 
faster than the Euler method because the integrand seems to be smoother. If one 
rejects our way to determine the necessary number of trajectories M as too crude 
and insists on taking the same number for both methods, the Ninomiya-Victoir 
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Method 


K 


M 


Rel. Error 


Time 


Euler 


32 


8192000 


0.00174 


91.94 sec 


Ninomiya-Victoir 


4 


2048000 


0.00204 


13.93 sec 


NV with drift 


4 


1024000 


0.00104 


2.88 sec 



Table 2. Computational time for the generalized SABR model 



method will require almost the same time as the Euler method in order to give 
comparable results at this level. However, the Ninomiya-Victoir method with drift 
still retains a convincing speed-up, again probably due to the simpler structure of 
the subroutines. 

4.3. Multi-dimensional generalized SABR. In this section we give the numeri- 
cal results corresponding to the multi-dimensional generalized SABR model. Recall 
that from Section 13.41 this model is given by 

dJQ(t) =a i Y« i X? i d§i 

AY.it) =Ki(6i - Yi)dt + hYidWi 

i = 1, . . . , N. For our experiment we choose N = 4 and as the derivative we take a 
basket option with the same weight on each stock. The parameters for the experi- 
ment have been chosen as follows: a = (1, 0.5, 0.3, 0.7) T , b = (0.5, 0.8, 0.4, 0.6) T , a = 
(0.5, 1, 0.7, 0.8) T , P = (0.6, 0.7, 0.8, 0.9) T , k = (0.2, 0.7, 0.5, 0.9) T , 6 = (0.3, 0.4, 0.6, 0.2) T , 
and, finally, 



P 



f 


1 


0.0111 


0.6395 


-0.1081 


-0.3414 


-0.0642 


-0.2054 


-0.0236\ 




0.0111 


1 


0.2698 


0.2770 


0.1651 


-0.3504 


-0.8186 


-0.4383 




0.6395 


0.2698 


1 


-0.1381 


-0.1379 


-0.0031 


-0.3169 


-0.0161 




-0.1081 


0.2770 


-0.1381 


1 


0.7312 


-0.9030 


0.0419 


-0.8121 




-0.3414 


0.1651 


-0.1379 


0.7312 


1 


-0.5969 


0.0747 


-0.6703 




-0.6420 


-0.3504 


-0.0031 


-0.9030 


-0.5969 


1 


0.1878 


0.8790 




-0.2054 


-0.8186 


-0.3169 


0.0419 


0.0747 


0.1878 


1 


0.2796 




-0.0236 


-0.4383 


-0.0161 


-0.8121 


-0.6703 


0.8790 


0.2796 


1 / 



Note that the above choice of p implies that B l and W l are negatively correlated, 
as usual in equity modeling. Moreover, p is positive definite - and in fact, chosen 
at random among all such correlation matrices. The estimated 'true value' of the 
basket option is 0.09254183. 

The convergence rates of the three different schemes for this experiment are 
graphically displayed in Figure |3J The picture is similar as in the two previous 
cases in the sense that there is second-order convergence for the two NV schemes 
and first-order convergence for the Euler scheme. 



Method 


K 


M 


Rel. Error 


Time 


Euler 


32 


2048000 


0.000934 


246.65 sec 


Ninomiya-Victoir 


4 


1024000 


0.002017 


52.33 sec 


NV with drift 


4 


1024000 


0.000862 


35.31 sec 



Table 3. Computational time for the multi-dimensional general- 
ized SABR model 



20 



CHRISTIAN BAYER, PETER FRIZ, RONNIE LOEFFEN 



o 



CD 



LU 



O 




CD 
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10 



20 



50 



100 



Number of timesteps 



Figure 3. Order of convergence for the multi-dimensional gener- 
alized SABR model. 

Further, the computational time for the generalized SABR model with a four- 
dimensional stock market, reported in Table |4~31 again shows the usual picture. The 
classical Ninomiya-Victoir method gives a speed-up between factors two and four 
(depending on the trust of the choice of M). In this case, one might, however, note 
that the error from the classical Ninomiya-Victoir method is more than twice higher 
the the errors from the competing methods. And again, the simpler structure of 
the ODEs in the case of a Ninomiya-Victoir method with drift leads to a convincing 
speed-up as compared to both other methods. 



where < /3 < 1, x >0 and h : [0, oo) — > M is continuous and either positive valued 
or negative valued. One can easily check that 



with a+ := max{a,0}, is a solution to ((8]). We briefly provide some details about 
uniqueness of the solution. When h takes only negative values, then the right hand 
side of the ODE ([5} is decreasing in the state variable and uniqueness follows for 
any x > (see e.g. Example 2.4 on p. 286 of [6]). When h takes only positive values 
and x > 0, uniqueness follows by the Picard-Lindeldf theorem. When h takes only 
positive values and x = 0, the solution ([9]) is not unique; for instance, y(t) = 
forms another solution. For this particular case, we have chosen, throughout the 



Appendix 



(8) 



In Sections [HU [2?2] and the following ODE appears: 

y'(t) =h{t){y{t)f, 
1/(0) =x, 



(9) 




t > 0, 
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paper, to work with the solution $o(t), since the flow $ x (t) is (right)-continuous 
at x = for all t > 0. 
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